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INTRODUCTION 


The  work  described  in  this  report  was  performed  in  collaboration  with  Professor 
Surendra  Singh,  who  came  to  the  Naval  Air  Warfare  Center  Weapons  Division 
(NAWCWD),  China  Lake,  California,  during  June  and  July  2005.  Professor  Singh,  a 
faculty  member  in  the  Electrical  Engineering  Department,  University  of  Tulsa,  Tulsa, 
Oklahoma,  was  visiting  as  an  Office  of  Naval  Research  -  American  Society  of  Electrical 
Engineers  (ONR-ASEE)  Summer  Faculty  Fellow.  He  worked  in  the  Optics,  RF,  and 
Material  Physics  Section  of  the  Physics  and  Computations  Sciences  Branch.  This  work 
involves  implementing  the  stabilized  version  of  the  bi-conjugate  gradient  algorithm  in 
order  to  iteratively  solve  a  linear  system  of  equations  resulting  from  integral  equations 
arising  in  electromagnetic  scattering  problems.  An  implementation  of  a  convergence 
acceleration  transform  to  speed  up  the  convergence  of  an  iterative  scheme  is  included 
here  as  the  appendix  to  the  report. 

The  solution  of  a  system  of  linear  equations  can  be  obtained  from  direct  methods, 
such  as  matrix  inversion,  and  indirect  methods  that  include  a  variety  of  iterative  schemes. 
The  conjugate  gradient  method  is  one  such  iterative  method  that  provides  an  additional 
advantage:  the  coefficient  matrix  need  not  be  stored  in  its  entirety  in  memory.  This 
provides  a  much  needed  flexibility  in  very  fine  discretizations  of  the  geometry  under 
investigation  without  taxing  the  memory  requirements.  We  implement  the  stabilized 
version  of  the  conjugate  gradient  method  and  provide  some  convergence  studies.  Noted 
that  the  stabilized  version  converges  very  rapidly  to  achieve  a  specific  precision,  and  then 
tends  to  oscillate  if  higher  precision  is  needed.  To  overcome  this  situation,  we  implement 
a  transform  to  accelerate  the  convergence  of  the  iterative  scheme.  The  results  of  applying 
this  transform  to  three  different  types  of  integral  equations  are  provided. 


BI-CONJUGATE  GRADIENT  STABILIZED  METHOD  (bcg-stab) 


The  bcg-stab  method  is  an  iterative  technique  to  solve  a  system  of  linear  equations  of 
the  form,  Ax  =  b .  The  matrix,  A ,  is  a  known  matrix  of  order  NxN ,  where  N  is  the 
number  of  unknowns,  b  is  the  known  forcing  function  of  length  N ,  and  x  is  the 
solution  vector  (length  N ).  The  method  can  be  implemented  without  storing  the  matrix 
A  .  The  algorithm  can  be  coded  such  that  only  a  row  or  a  column  of  the  matrix  is  needed 
at  a  time.  This  can  be  very  helpful  in  cases  where  storing  the  entire  matrix  would  pose  a 
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significant  storage  problem.  The  drawback  of  not  storing  the  matrix  is  that  it  gets 
computed  twice  in  each  iteration  of  the  method.  So  we  gain  advantage  in  storage  but  we 
pay  a  price  in  increased  computation  time. 

Here  is  the  pseudo-code  for  the  bcg-stab  method  (Reference  1): 

Initialization:  alpha=l;  rerr=100;  conv=10“4 ;  r0  =  b ;  u0  =  p0  =  0 ; p0  =  a0  =  co0  =  0  ; 

nitm^munkn; 
vall=dot(b,b) 
for  i  =  1 :  nitm 

if  (rerr  >  conv) 

P,  =  dot{b,r,_ ,) 

P  =  (P,  I Pt-i )(«/-■ 

P,  =  f,-\  +  P(P,-\  VM) 

v,  =  Apj  Note:  This  operation  can  be  split  up  so  that  we  use  one  row  at  a  time. 
at  =  p,  /  dot(b,  v, ) 
s  =  ri- 1  -  <*v, 

t  =  As  Note:  This  operation  can  be  split  up  so  that  we  use  one  row  at  a  time. 
coi  =  dot(t,s)/  dot(t,t) 
x,  =  xM  +  aipt  +  (o,s 
r=s-  (Ojt 

rerr  =  abs{dot{rj ,  rt )  /  vaJA) 

end 

end 


SUBROUTINE  bcgstab(b,nunkns, nitm, conv, x, nit, rerr) 

Input  Variables 

b  (complex  vector):  Right-hand  side  or  excitation  vector  of  length  nunkrts. 
nunkns  (integer):  Number  of  unknowns  (N), 

nitm  (integer):  Maximum  number  of  iterations  for  bcg-stab  (typically  set  equal  to 
mmlcns ). 

conv  (real):  Convergence  factor  to  stop  computation  for  bcg-stab  (typically  10  to  10  ). 

Output  Variables 

x  (complex):  Solution  vector  of  length  nunkns. 
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nit  (integer):  Number  of  iterations  taken  by  bcg-stab  method  to  converge. 

rerr  (real):  Residual  error.  The  bcg-stab  method  will  stop  when  the  residual  error  (rerr) 

becomes  less  than  the  convergence  factor  ( conv ). 

Lippman-Schwinger  Integral  Equation 

We  implement  the  bcg-stab  method  in  the  numerical  solution  of  the  Lippman- 
Schwinger  integral  equation  (Reference  2): 


E(r)  =  E  mc(r)  +  k20 


Jg  8  (r,  r')  •  Af(r')  E(r')  dr' 


(1) 


In  Equation  1,  the  symbols  have  their  usual  meaning:  r  =  (x,y) ,  E  is  the  total  electric 

field,  E""7  is  the  incident  electric  field,  Ga  is  the  two-dimensional  (2-D)  free-space 
Green  tensor  and  its  explicit  form  is  given  in  Reference  2,  k0=2jtlA,  and 
A£(r)  =  e  —  1 ,  when  r  is  entirely  within  the  nanowire  with  permittivity  e  ,  and  0 
otherwise.  The  nanowire  is  assumed  to  be  infinite  in  the  z  direction  and  all  spatial 
variations  occur  in  the  x-y  plane.  The  integration  is  carried  out  over  the  cross  section 
of  the  nanowire  that  is  embedded  in  vacuum.  We  solve  this  equation  numerically  for  a 
nanowire  30  run  in  size  (side  dimension)  illuminated  by  an  incident  field  of  wavelength, 
A  -  400  nm,  with  the  electric  field  vector  perpendicular  to  the  axis  of  the  nanowire.  In 
this  example  we  solved  for  zero  absorption,  that  is,  the  imaginary  part  of  the  permittivity 
is  zero. 

Table  1  convergence  comparison  of  beg  and  the  bcg-stab  for  two  values  of  the 
convergence  factor,  nrx  is  a  geometry  discretization  parameter  such  that  the  number  of 
unknowns  is  N  =  2* nrx* nrx.  #  is  the  number  of  iterations  taken  by  the  method  to 
converge. 


Table  1.  Convergence  Comparison,  nvx  and  N. 


Parameter 

No.  of 
unknowns 

conv  =  1 0  4 

conv  = 1 0  J 

nrx 

N 

beg  no. 

bcg-stab  no. 

beg  no 

bcg-stab  no. 

20 

800 

90 

57 

127 

174 

30 

1800 

320 

37 

424 

571 

40 

3200 

589 

19 

777 

688 
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Figures  1  and  2  show  the  convergence  behavior  of  the  bcg-stab  method  as  a  function 
of  the  iterations  for  conv  =  10  '4  and  conv  =  lO-5 ,  respectively.  We  see  that  the  method 
converges  extremely  rapidly  for  conv  =  10”4  in  comparison  to  beg,  which  takes  589 
iterations  to  converge,  as  indicated  in  Table  1.  The  numerical  results  shown  in  Table  1 
were  obtained  on  a  Pentium  IV  32-bit  PC.  The  performance  and  precision  of  the  bcg-stab 
method  may  improve  if  the  code  is  run  on  a  64-bit  machine. 
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FIGURE  1 .  Convergence  of  bcg-stab  Method  Showing  Log  of  Residual 
Versus  Iterations  for  conv  =  10‘4  and  nrx  =  40. 
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FIGURE  2.  Convergence  of  bcg-stab  Method  Showing  Log  of  Residual 
Versus  Iterations  for  conv  =  10  s  and  nrx  =  40. 


When  the  imaginary  part  of  permittivity  is  non-zero,  the  convergence  of  the 
algorithms  is  quite  rapid.  Table  2  is  convergence  comparison  beg  and  beg  (bcg__stab)  as  a 
function  of  wavelength,  X ,  and  convergence  factor  (conv).  The  results  are  obtained  from 
FORTRAN  code:  eigen2.  The  silver  nanowire  is  solid  and  the  the  imaginary  part  of 
epsilon  (ei)  &  0,  nrx  =  40,  (and  N  =  3200).  The  convergence  of  beg  stab  as  a  function 
of  the  iterations  is  shown  in  Figures  3  and  4. 


Tablel.  Convergence  Compairson,  A.  and  conv. 


X  —  300  run 
Number  of  iterations 

X  =  400  nm 

Number  of  iterations 

X  =  500  nm 

Number  of  iterations 

conv 

beg 

beg  stab 

beg 

bcg stab 

beg 

beg  stab 

10-6 

8 

5 

66 

41 

152 

83 

10~8 

11 

7 

101 

75 

249 

205 
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FIGURE  3.  Convergence  of  bcg  stab  Method  for  nrx  =  40  (Number  of 
Unknowns  =  3200),  X  =  400  run.  Convergence  Factor  =  10  8  (Result 
Obtained  From  MATLAB  Version  of  the  Code). 
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iteration  number 

FIGURE  4.  Convergence  of  bcg  stab  Method  of  nrx=  40  (Number  of 
Unknowns  =  3200),  X-  500  nm,  Convergence  Factor  (conv)  =  10‘8 
(Result  Obtained  From  MATLAB  Version  of  bcg  stab). 


LEVIN  TRANSFORM 


As  we  have  seen  in  the  results  shown  in  Figure  2,  the  bcg-stab  method  oscillates  for  a 
long  time  before  converging.  Our  goal  is  to  apply  a  transform  such  that  the  sequence  of 
vectors  given  by  the  iterative  scheme  (bcg-stab)  will  converge  faster  to  the  solution,  x ,  of 
the  system  Ax  -b .  There  are  a  number  of  convergence  acceleration  methods,  including 
Levin’s  transform,  Wynn’s  epsilon  algorithm,  Chebyschev-Toeplitz  algorithm,  and 
Brezinski’s  6  algorithm.  Here  we  focus  on  the  Levin  transform  that  provides  significant 
enhancement  in  convergence  of  a  vector  sequence.  To  begin  with,  we  start  with  the 
transform  as  applied  to  a  scalar  sequence  and  then  extend  it  to  be  applicable  to  a  vector 
sequence.  Let  Sn  be  the  partial  sum  of  n  terms  of  a  series  such  that  Sn  ->  S  as  n  -»  oo , 
where  S  is  the  sum  of  the  series.  If  the  partial  sums  are  coming  from  an  infinite  series 
that  converges  extremely  slowly,  then  the  scalar  sequence  S0,  S',,  S2 ,  ...,  will  take  a 
long  time  to  converge.  This  is  when  we  employ  techniques  to  enhance  the  rate  of 
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convergence  to  arrive  at  S  with  only  a  small  value  of  n .  To  this  end,  the  klh  order  Levin 
transform,  tk(Sn),  which  provides  an  estimate  of  the  sum  of  the  series,  S ,  may  be 
computed  as  (Reference  3): 


/(«)  - 
lk  “ 


n  +  i 


\n  +  k. 


V  ^n+i+l  _  ^n+i  J 


n  +  k  ) 


v  ^rt+i+i  Sn+I 


\ 


k  =  0,1,2,... 


(2) 


A  significant  advantage  of  the  Levin  transform  is  that  the  higher  order  iterates  are 
computed  from  the  partial  sums,  S0,  S, ,  S2 ,  rather  than  the  lower  iterates  as  in  the 

case  of  Wynn’s  algorithm,  thereby  providing  some  immunity  from  the  accumulation  of 
round-off  errors. 

Consider  the  solution  of  the  following: 


Ax  =  b  (3) 

where  A  is  a  known  matrix  (moment)  of  order  N  x  N  ;  A  -  [a.  ] ,  and  x  is  the  unknown 
column  vector  of  order  Axl;  x  =  [x,  x2,...xN]r ,  and  b  is  the  known  forcing  function 
column  vector  of  order  N  x  1 ;  b  =  [b]  b2 . .  ,bN ]7 .  With  A  and  b  known,  we  calculate  the 
initial  estimate  of  x ,  designated  as  jc  by  the  use  of  Gauss-Seidel  relaxation,  thus  arriving 
at  vectors  {S'q},  {*S'1 } ,  {S2  } ,  and  so  on.  Each  of  these  vectors  is  of  order  N  x  1 .  Once  we 
have  three  initial  vectors,  we  can  apply  the  Levin  transform  using  Equation  1  to  obtain  an 
estimate  of  the  solution  vector,  x  .  A  numerical  example  of  this  process  to  a  scalar 
sequence  is  given  in  Reference  4.  The  application  of  Equation  1  to  a  vector  sequence 
requires  a  bit  of  bookkeeping,  as  the  scalar  Levin  transform,  ,  now  becomes  a  vector, 

,  of  length  N  .  As  with  any  iterative  scheme,  the  computation  process  should  stop 
when  a  convergence  criterion  is  met.  We  employ  the  following  measure: 


||(^x  -  b) | 


<C, 


(4) 
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The  computations  are  stopped  when  the  above  criterion  is  met  for  a  pre-defined  convergence 
factor,  Cf ,  and  we  have  the  estimate  of  the  solution,  x  .  We  should  point  out  that  the  initial 

vectors  can  be  generated  from  a  variety  of  schemes.  It  may  be  possible  to  use  the  iterations 
from  a  beg  algorithm  and  then  subsequently  use  the  Levin  transform  or  any  other  vector 
acceleration  technique  to  enhance  the  convergence.  One  significant  advantage  of  using  the 
beg  estimates  is  that  the  moment  matrix,  A ,  does  not  need  to  be  stored  entirely  in  memory, 
as  the  algorithm  can  be  implemented  in  such  a  way  that  only  a  row  or  column  at  a  time  is 
needed  in  the  computations.  This  alleviates,  to  a  great  extent,  any  storage  problems 
routinely  encountered  when  solving  a  very  large  system  of  equations. 


NUMERICAL  EXAMPLES 


For  our  first  example,  we  consider  the  solution  of  the  integral  equation  for  the  charge 
distribution,  q(y),  on  a  conducting  strip  of  width  =  2  w  with  excitation,  V(y) : 


(5) 


For  illustration  purposes,  we  take  the  permittivity  of  the  medium,  £  = 1  F/m.  The 
integral  equation  is  cast  into  a  system  of  linear  equations  of  the  form  given  in  Equation  2 
by  using  method  of  moments  (MOM)  with  pulse  basis  functions  and  point  matching. 
Figure  5  shows  the  charge  distribution  for  N  =  150,  and  a  forcing  function,  V(y)  -  y2 . 
Note  that  the  Levin  transform  result  matches  very  well  with  direct  matrix  inversion. 
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FIGURE  5.  Charge  Distribution  on  a  Conducting  Strip  of  Half-Width  w  =  1 
m,  Forcing  Function  V(y)  =  y2 9  Convergence  Factor  C f  =  10“3,  Number  of 
Unknowns  N  =  150. 


The  second  example  to  illustrate  the  Levin  transform  involves  a  thin  wire  antenna. 
The  current  distribution  on  a  thin  wire  can  be  computed  by  the  numerical  solution  of 
Hallen’s  equation  or  Pocklington’s  equation.  Pocklington’s  integral  equation  for  the 
current  distribution,  I'(y)9  on  a  thin  conducting  wire  of  radius  =a,  length  =2 h,  oriented 
along  the  y  -  axis  is  given  by: 


ri 

\l(y')K{y-y')dy'  =  E‘y{y),  y  e  (~h,h) 


(6) 


where 
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-jk„R 


£(y-y')=  f  — d(t>’ 

2n  J  » 


(7) 


Here  the  distance  between  the  observation  and  source  point,  R ,  can  be  written  as: 


7?  =  [(T-/)2+4fl2sin2 


2 


(8) 


The  incident  electric  field  is  given  as  E'y (y)  =  VS(y-  yg)  for  the  antenna  case,  where 
V  is  the  feed  voltage  and  y  is  the  feed  point.  In  Equation  6,  rj  is  the  intrinsic 
impedance  and  kQ  is  the  wave  number  of  the  medium.  Once  again,  with  the  application 
of  the  method  of  moments,  with  suitable  basis  and  testing  functions,  we  can  cast  Equation 
6  in  the  form  of  a  system  of  linear  equations:  Ax-b.  The  current  distribution  on  a  center- 
fed,  half- wavelength  (h-/ 1/4)  thin  wire  antenna  with  V —  1  volt,  ^  =  0m, 

a  =  0.007022/1 ,  Cf  =10~3,  and  X  —  \  m  is  shown  in  Figure  6.  A  detailed  account  of 
these  examples  can  be  found  in  Reference  4. 
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FIGURE  6.  Current  Distribution  on  a  Center-Fed,  Half-Wavelength,  Thin 
Wire  Antenna  of  Radius  =  0.007022  X ,  Voltage  V  =  1  Volt,  Number  of 
Unknowns  N  =  21 ,  and  Convergence  Factor  Cf  =  10~3 . 
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Appendix 

FORTRAN  CODE  FOR  BI-CONJUGATE  GRADIENT  STABILIZED  METHOD 


FORTRAN  listing  of  SUBROUTINE  bcg-stab: 
subroutine  bcg_stab (b, nunkns, nitm, conv, ci, nit, rerr) 
implicit  none 

complex,  intent (in),  dimension ( 1 : nunkns)  ::  b 

complex,  intent (out),  dimension ( 1 : nunkns )  ::  ci 

integer,  intent (in)  ::  nunkns, nitm 

integer,  intent (out)  ::  nit 

real,  intent (in)  ::  conv 

real,  intent (out)  ::  rerr 

complex, dimension (1 : nunkns)  ::  p,r,v,s,t 
complex  ::  beta, alpha, vail, rho, rhoi, omega 
integer  : : i 
complex  : :  CDOTU 

ci=0 ; p=0 ; v=0 ; alpha=l ; omega=l ; rho=l ; r=b; s=0 ; t=0 
rerr=100 

vallr=dot_product  (b,  b) 
do  nit=l,nitm 

if  (rerr>conv)  then 
rhoi=dot_product (b, r) 

beta= (rhoi/rho) * (alpha/omega) 

p=r+beta* (p-omega*v) 

do  i=l, nunkns 

v (i) =CDOTU (nunkns, arow (i) , 1, p, 1) 
enddo 

alpha=rhoi/dot_product (b, v) 
s=r-alpha*v 

do  i=l, nunkns 
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t ( i ) =CDOTU ( nunkns , arow ( i ) , 1 , s , 1 ) 
enddo 

omega=dot_product ( t , s ) /dot_product ( t , t ) 

ci=ci+alpha*p+omega*s 

r=s-omega*t 

rho=rhoi 

rerr=abs (dot_product (r, r) /vail) 

!  print  * , ' iteration#= ' , nit , ' residual= 1 , rerr 
else 
exit 
end  if 
end  do 

end  subroutine  bcg_stab 
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